################################################################## 
## 全局变量

export work_dir=/public/home/xxf2019/20230331_colrNew/public_data/liver
export dna_preprocess_scripts=~/dna_preprocess_scripts
export config_path=${work_dir}/config
export config_file=${config_path}/config.sh

source ${config_path}/config.sh

################################################################## 
## 突变质控
## PASS
for file in `ls ${vcf_path} | grep filterStep1 | grep vcf.gz | grep -v tbi`
do
echo $file
tumor_normal=`echo ${file} | sed 's/.filterStep1.vcf.gz//'`
zcat ${vcf_path}/${file} | awk -F'\t' '{if($0~"#" || $7=="PASS")print}' > ${vcf_path}/${tumor_normal}_PASS.vcf
done

## Functator注释
cat ${config_path}/tumor_normal.list | grep -v Normal | awk -F, '{print $1,$2}' | xargs -P 10 -I sample sh -c '
sh ${dna_preprocess_scripts}/GATK4_Functator.sh sample ${config_path}
'

################################################################## 
## 新增质控
#### 3、提取需要质控的参数，见《20220225_体细胞突变位点质控流程.xlsx》
export cpu=20
for samples in ` cat ${config_path}/tumor_normal.list | grep -v Normal | tail -146` 
do
echo $samples
Tumor=` echo $samples | awk -F, '{print$1}' `
Normal=` echo $samples | awk -F, '{print$2}' `	
#sh ${dna_preprocess_scripts}/mutect2_Scripts/GATK4_SomaticVariant_Qc_Run.sh ${Tumor} ${Normal} ${config_file} ${cpu}
echo " sh ${dna_preprocess_scripts}/mutect2_Scripts/GATK4_SomaticVariant_Qc_Run.sh ${Tumor} ${Normal} ${config_file} ${cpu} " | \
qsub -N ${Tumor}_"SomaticQc" -l nodes=1:ppn=${cpu},mem=10gb,walltime=240:00:00 -q batch -d ${Qsub_log_path}
done

## 记录每个样本是否所有突变都在qc中有
echo "Tumor,Normal,Raw_vcf_muts,SomaticVariant_Qc_muts" > ${vcf_qc_path}/Record_SomaticVariant_Qc.csv

for samples in ` cat ${config_path}/tumor_normal.list | grep -v Normal | tr '\t' ',' ` 
do
echo $samples
Tumor=` echo $samples | awk -F, '{print$1}' `
Normal=` echo $samples | awk -F, '{print$2}' `
vcf_muts=`cat ${vcf_path}/${Tumor}_${Normal}_PASS.vcf | grep -v "#" | wc -l`
qc_muts=`cat ${vcf_qc_path}/${Tumor}_${Normal}_MutQc.csv | grep -v POS | wc -l`
echo "${Tumor},${Normal},${vcf_muts},${qc_muts}" >> ${vcf_qc_path}/Record_SomaticVariant_Qc.csv
done

####################################
## 检查数量一致性
## 若不一致，则重新计算不一致的突变
## 此步骤很重要，跑完再重跑上面计算突变数目
for Tumor in `cat ${vcf_qc_path}/Record_SomaticVariant_Qc.csv | awk -F',' '{del_num=($3-$4);if(del_num > 0){print $1}}'` 
do
Normal=` cat ${config_path}/tumor_normal.list | grep ${Tumor} | awk -F, '{print $2}'`

cat ${vcf_path}/${Tumor}_${Normal}_PASS.vcf | grep -v "#" | awk -F'\t' '{OFS=":"}{print $1,$2,$4,$5}' \
> ${vcf_qc_path}/tmp/${Tumor}_${Normal}_lossMut.check.Raw.tsv
cat ${vcf_qc_path}/${Tumor}_${Normal}_MutQc.csv | grep -v POS | awk -F',' '{print $1}' \
> ${vcf_qc_path}/tmp/${Tumor}_${Normal}_lossMut.check.mutqc.tsv
## 提取质控缺少的突变
diff_line=`cat ${vcf_qc_path}/tmp/${Tumor}_${Normal}_lossMut.check.Raw.tsv ${vcf_qc_path}/tmp/${Tumor}_${Normal}_lossMut.check.mutqc.tsv | \
sort | uniq -u`
diff_line_num=`cat ${vcf_qc_path}/tmp/${Tumor}_${Normal}_lossMut.check.Raw.tsv ${vcf_qc_path}/tmp/${Tumor}_${Normal}_lossMut.check.mutqc.tsv | \
sort | uniq -u | wc -l`
## 重新计算之前未计算的行
if [ ${diff_line_num} -gt 0 ]
then
for line in `echo ${diff_line} | tr ' ' '\n'`
do
echo $line
sh ${dna_preprocess_scripts}/mutect2_Scripts/GATK4_SomaticVariant_Qc_Step.sh ${Tumor} ${Normal} ${line} ${config_file} 
done
else
echo "The line is matched"
fi
##
done

####################################
#### 4、注释重复区域的信息
## 产生的文件如下
## ${vcf_qc_path}/${Tumor}_${Normal}_MutQc.annotationRegion.csv
for samples in ` cat ${config_path}/tumor_normal.list | grep -v Normal | tr '\t' ',' ` 
do
echo $samples
Tumor=` echo $samples | awk -F, '{print$1}' `
Normal=` echo $samples | awk -F, '{print$2}' `

sh ${dna_preprocess_scripts}/mutect2_Scripts/GATK4_SomaticVariant_Qc_AnnotationRegion.sh ${Tumor} ${Normal} ${config_file}

done

####################################
#### 5、整合MAF文件和对应的质控文件，完成质控
## 每一步质控位点的数目：${vcf_qc_path}/${Tumor}_${Normal}_MutQc.RecordQcNum.csv
## 满足质控要求的位点数目：${vcf_qc_path}/${Tumor}_${Normal}_MutQc.KeepMut.csv
for samples in ` cat ${config_path}/tumor_normal.list | grep -v Normal | tr '\t' ',' ` 
do
echo $samples
Tumor=` echo $samples | awk -F, '{print$1}' `
Normal=` echo $samples | awk -F, '{print$2}' `

sh ${dna_preprocess_scripts}/mutect2_Scripts/GATK4_SomaticVariant_Qc_Filter.sh ${Tumor} ${Normal} ${config_file}
done


## 记录所有样本每一步质控的reads
cat ${vcf_qc_path}/*MutQc.RecordQcNum.csv | grep -m 1 INDEL_Noncoding_TRF_Reatmasker| \
awk -F'\t' '{OFS=","}{print "Tumor,Normal,"$line }' > ${vcf_qc_path}/Record_All.MutQcNum.csv

for samples in ` cat ${config_path}/tumor_normal.list | grep -v Normal | tr '\t' ',' ` 
do
echo $samples
Tumor=` echo $samples | awk -F, '{print$1}' `
Normal=` echo $samples | awk -F, '{print$2}' `

cat ${vcf_qc_path}/${Tumor}_${Normal}_MutQc.RecordQcNum.csv | grep -v INDEL_Noncoding_TRF_Reatmasker | \
awk -F'\t' '{OFS=","}{print Tumor,Normal,$line }' Tumor=${Tumor} Normal=${Normal} >> ${vcf_qc_path}/Record_All.MutQcNum.csv
done

####################################
##### 7、提取满足质控的位点
for samples in ` cat ${config_path}/tumor_normal.list | grep -v Normal | tr '\t' ',' ` 
do
echo $samples
Tumor=` echo $samples | awk -F, '{print$1}' `
Normal=` echo $samples | awk -F, '{print$2}' `

sh ${dna_preprocess_scripts}/mutect2_Scripts/GATK4_SomaticVariant_Qc_AbstarctQcMut.sh ${Tumor} ${Normal} ${config_file}
done

## 记录bed提取前后的突变数数目
echo -e "Tumor,Normal,KeepMutbed_muts,QcVcf_muts" > ${vcf_qc_path}/Record_QcVcf.csv
for samples in ` cat ${config_path}/tumor_normal.list | grep -v Normal | tr '\t' ',' ` 
do
echo $samples
Tumor=` echo $samples | awk -F, '{print$1}' `
Normal=` echo $samples | awk -F, '{print$2}' `

KeepMutbed_muts=`cat ${vcf_qc_path}/${Tumor}_${Normal}_MutQc.KeepMut.csv | grep -v Hugo  | wc -l`
QcVcf_muts=`cat ${vcf_path}/${Tumor}_${Normal}_QC.vcf | grep -v "#" | wc -l`

echo -e "${Tumor},${Normal},${KeepMutbed_muts},${QcVcf_muts}" >> ${vcf_qc_path}/Record_QcVcf.csv
done


####################################
##### 8、满足质控的位点进行注释
for samples in ` cat ${config_path}/tumor_normal.list | grep -v Normal | tr '\t' ',' ` 
do
echo $samples
Tumor=` echo $samples | awk -F, '{print$1}' `
Normal=` echo $samples | awk -F, '{print$2}' `

echo " sh ${dna_preprocess_scripts}/GATK4_Functator_QC.sh ${Tumor} ${Normal} ${config_path} " | \
qsub -N ${Tumor}_"Funcotator" -l nodes=1:ppn=1,mem=10gb,walltime=240:00:00 -q batch -d ${Qsub_log_path}
done


####################################
#### 9、统计每个样本的SNV和INDEL情况
echo -e "Tumor\tNormal\tSNV_num\tIndel_num" > ${qc_path}/snv_indel_num.list
for samples in ` cat ${config_path}/tumor_normal.list | grep -v Normal | tr '\t' ',' ` 
do
echo $samples
Tumor=` echo $samples | awk -F, '{print$1}' `
Normal=` echo $samples | awk -F, '{print$2}' `

snv_num=`cat ${maf_path}/${Tumor}_${Normal}_QC.maf | grep -w SNP | wc -l`
indel_num=`cat ${maf_path}/${Tumor}_${Normal}_QC.maf | grep -E -w "INS|DEL" | wc -l`

echo -e "${Tumor}\t${Normal}\t${snv_num}\t${indel_num}" >> ${qc_path}/snv_indel_num.list
done


####################################
#### 10、用到的突变数量
echo "Tumor,Normal,MutQcNum" > ${qc_path}/Vcf_QC.list
for samples in ` cat ${config_path}/tumor_normal.list | grep -v Normal | tr '\t' ',' ` 
do
echo $samples
Tumor=` echo $samples | awk -F, '{print$1}' `
Normal=` echo $samples | awk -F, '{print$2}' `
QcVcf_muts=`cat ${vcf_path}/${Tumor}_${Normal}_QC.vcf | grep -v "#" | wc -l`
echo ${Tumor}","${Normal}","${QcVcf_muts} >> ${qc_path}/Vcf_QC.list
done
